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ABSTRACT 

We study the dynamical interaction in which the two single runaway stars AE Aurigae 
and /i Columbae and the binary i Orionis acquired their unusually high space velocity. 
The two single runaways move in almost opposite directions with a velocity greater 
than 100 km s -1 away from the Trapezium cluster. The star l Orionis is an eccentric 
(e ~ 0.8) binary moving with a velocity of about lOkms -1 at almost right angles 
with respect to the two single stars. The kinematic properties of the system suggest 
that a strong dynamical encounter occurred in the Trapezium cluster about 2.5 Myr 
ago. Curiously enough, the two binary components have similar spectral type but very 
different masses, indicating that their ages must be quite different. This observation 
leads to the hypothesis that an exchange interaction occurred in which an older star 
was swapped into the original l Orionis binary. We test this hypothesis by a combina- 
tion of numerical and theoretical techniques, using N-body simulations to constrain 
the dynamical encounter, binary evolution calculations to constrain the high orbital 
eccentricity of i Orionis and stellar evolution calculations to constrain the age dis- 
crepancy of the two binary components. We find that an encounter between two low 
eccentricity (0.4<e£0.6) binaries with comparable binding energy, leading to an ex- 
change and the ionization of the wider binary, provides a reasonable solution to this 
problem. 
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1 INTRODUCTION 

OB runaways are a subgroup of spectral type O and B 
stars which are recognized by their unusually high velocity 
(>30kms _1 ), often away from known star forming regions 
(Blaauw, 1961). Stone (1979) also included stars at large dis- 
tance (up to several kpc) from the Galactic plane for which 
no parent association could be found. About 40 per cent 
of O stars and 5-10 per cent of B stars are runaways (Stone 
1991), and only one in ten OB runaways has a known binary 
companion (Gies & Bolton 1986). 

The peculiar velocities of about 50 OB runaways have 
been determined and the trajectories of ~ 20 can be traced 
back to a nearby association (Hoogerwerf et al. 2000, 2001), 
supporting the idea that these stars originated in stellar clus- 
ters from which they were later ejected with high velocities. 
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The two main mechanisms that have been proposed to 
explain the high velocities of OB runaways are: 

1. ejection upon a supernova in a binary system (Zwicky 
1957; Blaauw 1961); 

2. ejection by a close multiple encounter in a star cluster 
(Poveda, Ruiz & Allen, 1967). 

We refer to scenario 1 and 2 as supernova ejection and dy- 
namical ejection, respectively. 

In the supernova ejection scenario the OB runaway gets 
a high velocity from the explosion of the companion star. If 
the binary is dissociated by the supernova, the companion 
star is ejected with a velocity of the order of its orbital ve- 
locity while if the binary remains bound, the mass loss in 
the supernova induces a recoil velocity in the center of mass 
of the binary (Blaauw 1961; Hills 1980; Tauris & Takens 
1998). This scenario predicts that about 30 per cent of the 
OB runaways should still be in a binary with a compact 
companion. Large natal kicks (up to 1000 kms -1 ) are ob- 
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served in isolated pulsars (Lyne & Lorimer 1994; Cordes 
& Chernoff 1998), but are not sufficient to explain the low 
binarity fraction observed for OB runaways. Even if the for- 
mation of a neutron star is accompanied by a kick drawn 
from a Maxwellian distribution with a mean of 600 km s -1 , 
about 30 per cent of the binaries remain bound (Portegies 
Zwart 2000). Selection effects may be important in prevent- 
ing the observation of compact companions to the runaways 
(e.g. Colin et al. 1996), but the supernova ejection scenario 
seems insufficient to explain all the known OB runaways 
(Portegies Zwart 2000). 

The supernova ejection scenario, however, satisfactorily 
explains the high velocities (~50kms _1 ) of some known 
high- mass X-ray binaries (van den Heuvel et al. 2000). In one 
of these, Vela X-l, the high space velocity is inferred from 
the morphology of a bow shock (Kaper et al. 1997) while the 
neutron star is visible as an X-ray pulsar. The only known 
example of a binary dissociated by a natal kick is given by 
the 09.5V star < Oph and the pulsar PSR J1932+1059. 
The large rotational velocity (Penny 1996) and helium abun- 
dance (Herrero et al. 1992) of £ Oph indicate that a phase of 
mass transfer occurred, as expected for runaways ejected by 
a supernova explosion. Under this hypothesis, a natal kick of 
about 400 km s -1 can be inferred for the pulsar (Hoogerwerf, 
de Bruijne & de Zeeuw 2001). 

An alternative scenario to the supernova ejection is the 
dynamical ejection scenario (Poveda, Ruiz & Allen 1967; 
Gies & Bolton 1986; Leonard 1990), in which close encoun- 
ters between binaries and/or single stars in dense stellar 
regions can lead to the ejection of stars with high velocity. 
Young stellar clusters and OB associations in early evolu- 
tionary stages can have high stellar densities ( ~ 10 J pc~ 3 ) 
in their cores and therefore close encounters between single 
stars, binaries or higher order systems can be very frequent. 
These encounters can result in the ejection of stars with high 
velocities. O and B stars are more likely to participate in 
strong encounters than low mass stars because of mass seg- 
regation and gravitational focusing. As a consequence, the 
number of close encounters involving O and B type stars 
can be high despite their short lifetime, especially when the 
cluster is young and the binary fraction is high. 

To help understanding the runaway mechanism we de- 
fine the kinematical age rki n as the time since the ejection 
from the parent association. This time can be derived inte- 
grating the trajectory of the runaway in the Galactic poten- 
tial. At first order Tkin = d/v, where d is the distance trav- 
eled from the cluster and v the velocity of the runaway. The 
kinematical age can help to identify the ejection mechanism. 
A runaway which is ejected by a supernova explosion has a 
kinematical age smaller than the age of the parent cluster 
because the primary star in the original binary evolved for a 
few million years before it exploded ejecting its companion. 
Dynamical ejection is likely to occur when the cluster is very 
young and therefore in this case the kinematical age is often 
similar to the age of the parent system. 

An example of stars ejected by a dynamical encounter is 
provided by the single runaways AE Auriga and fi Columba 
and the spectroscopic binary t Orionis, whose trajectories in 
the Galactic potential have been traced back in time and 
intersect about 2.5 Myr ago at the expected location of the 
Trapezium cluster (Hoogerwerf et al. 2001). The runaways 
AE Auriga and fi Columbae show no evidence of binary evo- 



lution, like a high rotational velocity or an increased helium 
abundance, and move in almost opposite directions away 
from the Trapezium cluster. The kinematical and evolution- 
ary properties of the system suggest the hypothesis that the 
stars were involved in a 4-body encounter that ejected them 
with high velocity. 

In this paper we describe a combined use of direct 
N-body simulations and stellar and binary evolution cal- 
culations to identify the type of encounter that ejected 
AE Auriga, /i Columbae and t Orionis. The paper is orga- 
nized as follows: In section § 2 we describe the main obser- 
vational properties of the stars and a possible evolutionary 
model for i Orionis, AE Auriga and pi Columba and we pro- 
pose a binary-binary encounter able to reproduce the cur- 
rent positions, velocities and physical properties of the four 
stars; in section § 3 we describe the code used to perform 
simulations of 4-body encounters and the classification of 
the outcomes; in section § 4 we describe the choice of initial 
conditions for the simulations, we present the main results 
obtained from the scattering experiments and we compare 
them with the observed properties of the stars while in sec- 
tion § 5 we present aspects of the evolution of the l Orionis 
binary after the encounter that can affect the model for the 
dynamical encounter. 



2 THE RUNAWAY STARS AE AURIGAE, fx 
COLUMBAE AND THE RUNAWAY BINARY 
l ORIONIS 

2.1 The kinematics 

The single runaway stars AE Aur and yi Col are O type 
stars moving away from the Orion star forming region. 
Blaauw & Morgan (1954) recognized them as a remarkable 
pair of runaways having similar spectral types and mov- 
ing in almost opposite directions with space velocities of 
about 100 km s" 1 . At a distance of about 250 pc from the 
Trapezium cluster the kinematical age of both stars is about 
2.5 Myr. 

The same kinematical age and parent cluster for the two 
stars together with their relative straight angle of motion 
suggest a common origin in a dynamical encounter. Gies & 
Bolton (1986) advanced the hypothesis of a binary-binary 
interaction that ejected AE Auriga and p, Columba and left 
i Orionis as the surviving binary. They noticed that t Ori is 
one of the most massive objects in the Orion nebula, has a 
high eccentric orbit and its component stars have relative 
orbital velocities similar to the space velocities of AE Aur 
and jj, Col. 

Parallaxes and proper motions provided by Hipparcos 
together with radial velocity measurements (Turon et al. 
1992) are shown in Table 1 and allow precise determina- 
tions of the Galactic trajectories of the binary and the single 
stars. The integration of the trajectories of AE Aur, p, Col 
and i Ori (Hoogerwerf et al. 2001) in the Galactic potential 
shows that 2.5 ± 0.2 Myr ago all the stars had a minimum 
separation of about 4pc. The position and velocity of the 
parent cluster, derived from the orbit integration and the 
4-body center of mass velocity, are consistent with those of 
the Trapezium cluster in the Orion nebula. 

We here summarize the arguments in favor of a common 
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Table 1. Parallax, distance from the sun and proper motions 
from the Hipparcos Catalogue (ESA 1997), radial velocities from 
the Hipparcos Input Catalogue (Turon ct al. 1992), velocity rel- 
ative to the Local Standard of Rest and velocity relative to the 
center of mass of the 4-body system for AE Aur, /i Col and i Ori. 



Table 2. Properties for AE Aurigae, /iColumbae (Gies and Lam- 
bert 1992; Bagnuolo, Riddle, Gies 2001) and the primary (t Ori A) 
and secondary (t Ori B) component of i Orionis (Stickland ct al. 
1987; Bagnuolo, Gies and Penny 1994). 



AE Aur 



7r [mas] 
d [pc] 

Mc* [mas yr" 1 ] 
M [mas yr _1 ] 
V rad [kins" 1 ] 
Vlsr [kms -1 ] 
Vcm [kms -1 ] 



2.24 ± 0.74 
446±™ 

-4.05 ± 0.66 

43.22 ± 0.44 
57.5 ± 1.2 

111.9 ± 24.2 
115 ± 5 



H Col 

2.52 ± 0.55 
397j£° 

3.01 ± 0.52 
-22.62 ± 0.50 

109.0 ± 2.5 

108.2 ± 4.2 
103 ± 2 



i Ori 
2.27 ± 0.65 

2.27 ± 0.65 
-0.62 ± 0.47 
28.7 ± 1.1 
17.4 ± 0.8 
18 ± 1 





AE Aur 


H Col 


l Ori A 


t Ori B 


Spectral type 


09.5 V 


O9.5V/B0 


09 III 


B0.8 III-IV 


T cf f [°K] 


31420 


31790 


32500 


27000 


log(g) 


4.07 


3.85 


3.73 


3.78 


Mass [Mq] 


16 


16 


39 


19 


Radius [Rq] 


6.0 


8.0 


16 


11 


Period [days] 








29.134 


Eccentricity 








0.764 



origin of the runaways AE Aur, /x Col and the i Ori binary 
in the Trapezium cluster: 

1. The trajectories of the single stars and the binary in- 
tersect 2.5 Myr ago in the Trapezium cluster. 

2. The velocity of the center of mass of the 4-body system 
is compatible with the velocity of the Trapezium cluster. 

3. The ages of all the four stars (see section 2.2 and 5.2) are 
consistent with the range in ages of the stars in the cluster 
(Palla & Stahler 1999). 

4. The high stellar density (> 2 x 10 4 pc~ 3 , McCaughrean 
& Stauffer 1994) in the core of the Trapezium favors dynam- 
ical interactions. 



2.2 The stellar evolution 

The t Orionis binary consists of a 09 III primary and a B0. 8 
III-IV secondary with a mass ratio of q=0.5 (Stickland et al. 
1987) or g=0.57 (Marchenko et al. 2000). Using the single- 
star evolutionary tracks by Schaller et al. (1992), Bagnuolo 
et al. (2001) estimate a difference in age of roughly a factor 
of two in the components and conclude that the system did 
not co-evolve. Therefore they propose that the system orig- 
inated in a binary-binary encounter in which an exchange 
interaction occurred. The age difference cannot be explained 
by a phase of mass transfer because of the high eccentricity 
of the binary: mass transfer would have circularized the or- 
bit of the binary. Exchange encounters are known to alter 
the eccentricity of the interacting binaries and could pro- 
vide a natural explanation for the high eccentricity of the 
system. Observational data on i Ori, AE Aur and jj, Col are 
presented in Table 2. 

Mason et al (1998) list i Ori with a speckle companion 
at a separation of about 0.1 ". This could mean that i Ori 
is the primary component of a hierarchical triple system. If 
the speckle companion is indeed bound to the binary, it has 
considerable consequences for the proposed origin of i Ori, 
as we discuss in § 5.1. However, at the moment it is not clear 
whether the third component is bound to the binary. 

A possible evolutionary scheme for the system is shown 
in Table 3. Starting from the currently observed parameters 
reported in Table 2, the stellar and binary evolution package 
SeBa (see Portegies Zwart & Verbunt 1996) is used to infer 



Table 3. Scheme of the evolutionary model adopted in this work 
for i Ori, AE Aur and fi Col. Current values of mass, radius and 
age for AE Aur, fj, Col (Gies and Lambert 1992; Bagnuolo, Riddle, 
Gies 2001) and the primary (t Ori A) and secondary (t Ori B) 
component of i Ori (Bagnuolo, Gies and Penny 1994) are used to 
predict masses and radii at the moment of the encounter. These 
parameters are used as initial conditions in the simulations. 
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=0 


(i, Ori B, n Col) 








l Ori B 


22 Mq 


6.0 Rq 






At Col 


18 Mq 


5.0 Rq 


T= 


=3.5 Myr 


(t Ori A, AE Aur) 








i Ori A 


42 Mq 


8.5 Rq 






AE Aur 


18 Mq 


5.0 Rq 


T= 


=4.5 Myr 


(t Ori A, 


AE Aur) 4 


- (t Ori B, n Col) 






-> (i Ori 


A, i Ori B] 


) + AE Aur + p, Col 






i Ori A 


42 Mq 


10 Rq 






i Ori B 


21.5 Mq 


9.0 Rq 






AE Aur 


18 Mq 


5.5 Rq 






ii Col 


18 Mq 


7.0 Rq 


T= 


=7 Myr 


(i Ori A, 


i Ori B) + 


AE Aur + \i Col 






l Ori A 


39 Mq 


20 Rq 






l Ori B 


20 Mq 


16 Rq 






AE Aur 


18 Mq 


6.5 Rq 






H Col 


18 Mq 


10 Rq 



masses and radii of each single star at the moment of the 
encounter. 

2.3 The encounter 

We explore the hypothesis of a dynamical ejection for 
AE Aur, n Col and i Ori using numerical simulations of 
binary-binary encounters and try to infer the interaction 
that is needed to produce the observed properties of the 
stars. First of all, one binary must be ionized during the 
encounter in order to eject the two single stars. In addition, 
an exchange is needed to solve the age discrepancy in the 
i Orionis binary. An example of a dynamical encounter pro- 
ducing a system like i Ori, AE Aur and fi Col is shown in 
Figure 1. A binary consisting of i Ori B and fx Col inter- 
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Figure 1. Example of a binary-binary encounter leading to the 
ionization of a binary and the exchange of one star. The two 
binaries are initially in the top right and bottom left corner and 
approach in the direction indicated by the arrows. The binary 
containing fi Col is unbound and i Ori B is exchanged in the other 
binary. As a result, the tOrionis binary is formed and AE Auriga; 
and /^Columbae are ejected as singles. 




acts with a binary containing i Ori A and AE Aur. In the 
encounter the softer binary is unbound releasing li Col and 
l Ori B, which is exchanged in the other binary and together 
with i Ori A forms the i Orionis binary while AE Aur and 
fj, Col escape in almost opposite directions. 



3 BINARY-BINARY SCATTERINGS 

3.1 The code and the setup of initial conditions 

The numerical simulations described in this pa- 
per are carried out with the scatter package in- 
cluded in the STARLAB software environment 
(McMillan & Hut 1996; Portegies Zwart et al. 2001; 
http://www.ids.ias.edu/~starlab). An N-body en- 
counter is resolved integrating the equations of motions 
of all the particles under the influence of their mutual 
gravitational forces using a fourth-order Hermite predictor- 
corrector scheme (Makino & Aarseth 1992). 

We consider a target binary composed of stars of mass 
M t i and M t2 , semi-major axis a t and eccentricity e t , and a 
projectile binary composed of stars of mass M p \ and M p2 , 
semi-major axis a p and eccentricity e p . Standard N-body 
units (Heggie & Mathieu, 1986) are used in which M t i + 
Mt2=l, a t =l and G=l. The relative velocity Voo is given in 
units of the critical velocity V c for which the total energy of 
the system in the 4-body center of mass reference frame is 
zero: 

JO (MuM^ + M p1 M p2 \ 
p\ at a p J 



where fj, = (M t i + M t2 )(M pl + M p2 )/M tot is the reduced 
mass and M to t the total mass of the system. Additional pa- 
rameters like the orbital phases and the relative orientation 
of the two binaries are randomly drawn from uniform distri- 
butions (Hut & Bahcall 1983). The initial eccentricities are 
drawn from a thermal distribution P(e) = 2e (Heggie 1975). 

The impact parameter b is chosen between zero and 
a maximum value according to an equal probability distri- 
bution for b 2 . The maximum value 6 max is determined au- 
tomatically for each experiment. The code defines different 
cylindrical shells with the same cross-sectional area and per- 
forms an arbitrary number of scatterings in successive shells 
until all the encounters in the outermost shell result in a 
preservation of the binaries. The limiting impact parameter 
of the outermost shell defines the value of b max specific to 
the experiment under consideration. 

Energy conservation in our simulations is typically bet- 
ter than one part in 10 6 . The error in energy conservation 
is checked in each experiment and if it exceeds 10 -5 the en- 
counter is rejected. The accuracy in the code is chosen in 
such a way that at most a few per cent of the encounters 
are rejected. 

3.2 Classification system for binary-binary 
encounters 

When the integration of the equations of motion is stopped 
(see appendix B for the stopping criteria) the encounters are 
classified. During the calculation we keep track of the near- 
est neighbors of all the stars and we compute the distances 
between them. If two stars approach each other within a 
distance smaller than the sum of their radii we classify the 
encounter as a collision. 

Many outcomes are possible in binary-binary interac- 
tions: two bound systems can be left, or one binary and two 
single stars, or a triple and a single star, or four single stars. 
The encounters are classified as follows: (i) preservation or 
flyby if the two original binaries remain bound, even though 
strongly perturbed; (ii) exchange if one star is exchanged 
by the binaries and two new binaries are formed; (iii) ion- 
ization if one binary is unbound leaving a binary and two 
single stars recoiling to infinity; (iv) exchange-ionization if 
one binary is unbound and one star is exchanged in the other 
binary; (v) triple if one binary is unbound and one of the 
stars is captured by the other binary in a bound stable orbit 
while the other star is ejected as single; (vi) total ionization 
if both binaries are unbound and four single stars recede to 
infinity. 

In Figure 2 we present different possible outcomes (and 
corresponding Feynman diagrams) of binary-binary interac- 
tions. 



3.3 Testing the code 

In order to test our scattering package, we performed a se- 
ries of binary-binary scattering experiments and compared 
the results with those of Mikkola (1983). He used a hybrid 
code containing three different integration methods: the hi- 
erarchy method, where one or more of the relative motions 
in the system are integrated as a perturbed keplerian orbit, 
direct integration and Heggie's general N-body regulariza- 
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Figure 2. Examples of different outcomes in binary-binary scattering events and the corresponding Fcynman diagrams. In the diagrams 
the paths of two stars are placed close together when they are bound, and the wavy lines represent gravitational perturbations. 
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tion method (Heggie 1974). Mikkola considered two identi- 
cal binaries with equal mass stars approaching each other 
with various relative velocities. We classify his results in the 
following way: unresolved for an encounter that cannot be 
classified as any of the other types at the moment the inte- 
gration is stopped; flyby if two binaries are in a hyperbolic 
relative orbit; stable triple for a hierarchical triple and one 
escaping star; single ionization if one binary and two escap- 
ing stars remain and finally total ionization if both binaries 
are disrupted. The frequencies of the different outcomes are 
reported in Table 4 as a function of the relative kinetic en- 
ergy Too between the binaries together with the results by 
Mikkola. The agreement is in most cases better than the 1-cr 
level and always within 2-er. 



We also computed the cumulative cross-section for close 
approach during the encounter between two identical circu- 
lar equal-mass binaries and compared them with those of 
Bacon, Sigurdsson & Davies (1996) and Cheung et al. (in 
preparation). The comparison with the most recent results 
by Bacon et al. and with the results by Cheung et al. in- 
dicates that the cross-sections are consistent within the nu- 
merical uncertainties. 



To test our 4-body scattering package against the 3- 
body package we reproduced the results by Rasio, McMillan 
and Hut (1995) regarding the formation of the triple system 
PSR B1620-26 in M4. Branching ratios and cross-sections 
obtained simulating binary-binary encounters are well con- 
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Table 4. Comparison of the frequencies of different types of outcomes with the results by Mikkola (1983). For the scattering experiments 
we consider binaries with four equal-mass stars and equal semi-major axes. 
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sistent with the ones obtained in the case binary-single star 
encounters. 



Figure 4. Distribution of velocities relative to the center of mass 
of the 4-body system for t Ori, AE Aur and fi Col (solid line) 
compared with the observed values (dashed line) for a=2. 



4 NUMERICAL SIMULATIONS 
4.1 Initial conditions 

We consider interactions between a target binary composed 
by AE Aur and the primary component of i Ori (t Ori A) 
and a projectile binary composed by n Col and the sec- 
ondary component of t Ori (i Ori B). The values of masses, 
radii and ages used as initial conditions for the simulations 
are reported in table 3. After the two binaries and their rela- 
tive orbits are selected we integrate the equations of motion 
until the encounter is resolved (see appendix B). The rela- 
tive velocity at infinity between the centers of mass of the 
two binaries is set equal to the mean dispersion velocity in 
the Trapezium cluster (2kms _1 , Herbig & Tendrup 1986). 
Conservation of total energy and momentum in the 4-body 
system center of mass frame is applied to derive the bind- 
ing energies of the initial binaries. Only the sum of the two 
binding energies can be derived from the conservation laws, 
so the ratio between the two remains as a free parameter: 
a = E t /E p . 

For each value of a between 0.1 and 10 we perform 
2000 scattering experiments and the resulting branching ra- 
tios and cross-sections are presented in Figure 3. Most of 
the encounters lead to the disruption of one binary, possi- 
bly accompanied by an exchange or a capture in a stable 
triple. Only about 10 per cent of the encounters result in an 
exchange. The fraction of flybys decreases as a function of 
a while the fraction of triples increases for a > 1. A value 
a > 1 means that the projectile binary is softer than the 
target binary and is more easily ionized. The larger is a the 
wider is the projectile binary. After an ionization, one of 
the component stars can be captured by the target binary 
which contains the most massive star of the system. The 
normalized cross section decreases steadily as a function of 
a because of the decrease in the normalization factor. For 

values a>4 the cross-section scales as a~ 2 . 

We mainly focus on encounters of the type 

(t Ori A, AE Aur) + (t Ori B, fi Col) -» (i Ori A, i Ori B) + AE Aur + fj Col 

leading to the ionization of one binary and the exchange of 
one specific star. This interaction is most favored for a < 3 
(see Figure 3). For clarity we select a=2 and for this value 
we perform 15000 scattering experiments. 
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4.2 The point particle limit 

We first analyse the properties of the stars resulting from the 
simulations in the approximate case of point particles. The 
final velocity distributions for i Ori, AE Aur and /i Col are 
shown in Figure 4. The velocities are relative to the center 
of mass of the 4-body system and the observed values (see 
Table 1) are presented as dashed lines. Due to momentum 
conservation the two single stars recoil with a velocity ~ 3 
times higher than the binary. Typical velocities for the stars 
range from 30 to 100 km s -1 , with an average of the order 
of their orbital velocities (~60 — 70kms~ 1 ) in the initial 
binaries. The velocity of the binary has a peak at about 
25kms~ 1 , close to the observed velocity 18 ± lkms -1 . 

We observe a correlation between the velocity of the bi- 
nary and that of the single stars, though with a large scatter, 
while there is no apparent correlation between the space ve- 
locities of the two single stars (see Figure 5). This might 
be due to the many degrees of freedom present in 4-body 
scatterings. 

In Figure 6 we show the distribution of the semi-major 
axis of the binary after the encounter. The observed value 
(Stickland 1987; Marchenko et al. 2000) is consistent with 
the most probable value in the distribution. 

The distribution of velocities and semi-major axis does 
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Figure 3. Left: branching ratios for different types of outcome as a function of the ratio a between the initial binding energies of the two 
binaries. Only flybys, exchanges, ionizations, triples and encounters leading to systems like AE Aur, fi Col and i Ori are shown. Right: 
normalized total cross-section as a function of the parameter a. 





not change for a g 3 but become significantly different for 
larger values of a. The consistency of the distributions with 
the observed values reflects the right choice of the total en- 
ergy in the encounter. The total energy of the 4-body system 
is known at present and its conservation before and after 
the encounter was used in the choice of the initial condi- 
tions. Different choices of the total energy available in the 
interaction lead to different distributions for the velocities 
and the binary semi-major axis. This seems to exclude the 
possibility that a fifth body was involved in the encounter. 

In order to verify the possibility of producing two sin- 
gle runaways moving in almost opposite directions, as is ob- 
served for AE Aur and n Col, we study the distribution of 
the angular displacement of the two stars with respect to the 
center of mass of the 4-body system. We compute the angle 



9 between the velocity vectors of the two single stars and 
report its distribution fg in Figure 7. We find that about 35 
% of all encounters result in an angular displacement 9 in 
the 45 degrees range 135° < 9 < 180°. 

The eccentricity of the binary after the encounter is 
consistent with a thermal distribution and therefore more 
than 50 per cent of the encounters result in a binary with 
e > 0.67. 

We computed the absolute fractions and the cross- 
sections of different types of encounters using the procedure 
described in appendix A. If ax is the cross-section for pro- 
cess X, a normalized cross-section can be defined as 



ox = 
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Figure 6. Distribution of the final binary semi-major axis (solid 
line) compared with the observed value (dashed line) for a=2. 




0.4 0.6 
a (AU) 



Figure 7. Distribution of the relative angle 9 between the velocity 
vectors of AE Aur and fi Col (solid line) with respect to the center 
of mass of the 4-body system. The dashed line indicates the value 
of the relative angle derived from observations. 




180 
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where «oo represents the relative velocity between the cen- 
ters of mass of the interacting binaries. As shown in Table 
5, about 1 in 6 encounters results in the ionization of one 
binary and the exchange of one specific star, leading to sys- 
tems like AE Aur, fi Col and l Ori. 

From the cross-section it is possible to derive the typical 
time-scale of the process under consideration if the stellar 
density and the stellar velocity dispersion of the region are 
known 

i i 

(3) 



7t (a 2 + a%) V c 2 n f b a x ' 



where n represents the mean stellar density in the cluster 
and fb is the fraction of binaries in the core. Substituting 



Table 5. Branching ratios fx and normalized cross-sections ax 
for different types of outcomes in the case a=2. la error bars are 
of the order of 1 per cent. We specifically split off encounters which 
lead to a system like i Ori, at the bottom, below the horizontal 
line. 



Encounter 


fx 


ax 


Preservation 


0.21 


1.20 


Exchange 


0.03 


0.18 


Ionization 


0.17 


1.00 


Exchange- ionization 


0.35 


2.06 


Triple 


0.24 


1.42 


Total ionization 


0.00 


0.00 


i Orionis 


0.20 


1.17 



typical parameters of the Trapezium cluster (with a binary 
fraction fb ~ 0.6; Duquennoy & Mayor 1991) we obtain 



3Myr / 5 x 10 4 pc" 



<Jx 



AU 2 



v 

aoT 



5 kms 



where v = Voo/V c . In the case a — 2 the typical time-scale 
for a binary-binary encounter resulting in an exchange- 
ionization is about 3 Myr. 



4.3 Collisions in binary-binary encounters 

In this section we relax the hypothesis of point mass stars 
and investigate the effects of physical radii on the outcomes 
of binary-binary encounters. In the simulations described in 
the previous sections, we keep track of the minimum dis- 
tances between all stars. In this way we are able to use the 
same data and test the effect of finite stellar radii, assuming 
that stars collide if they approach each other within some 
minimum distance rf co ii- 

Since the initial eccentricities of both binaries are ran- 
domly drawn from a thermal distribution between and 
1, stars in very eccentric orbits may collide at the first peri- 
center passage. This would not be realistic, and therefore we 
limit the range in initial eccentricities. In a highly eccentric 
orbit two stars with total radius R collide if R > 0.75a(l — e) 
(Kopal 1959). Inverting this equation leads to our defini- 
tion of e max . For the observed stellar radii (see Table 3) 
e m ax = 0.9. We therefore exclude binaries with initial eccen- 
tricity e > 0.9 from further consideration. 

Figure 8 shows the cumulative distribution for mini- 
mum distances between any two stars for the ~ 8000 en- 
counters which remained after selecting the systems with 
e < e max (solid curve) and for circular initial binaries 
(dashed curve). Since circular orbits are very rare in our 
standard initial conditions, we performed a separate simu- 
lation of 8000 encounters. If we assume stellar radii from 
Tab. 3 at the moment of the encounter, a collision occurs if 
Tmin/a S3 0.07. In that case about half of the initially circular 
binaries result in a collision, whereas only 25 per cent of the 
eccentric systems survive merging. In figure 9 we show the 
fraction of binaries that experience a collision as a function 
of e max . The two figures 8 and 9 show similar trends; high 
initial eccentricities enhance the collision rate. 

In Figure 10 we show the proportion of ionized systems 
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Figure 8. Cumulative distribution for minimum distances r m ; n , 
between stars during binary-binary encounters. The solid line rep- 
resents the result for systems with e < 0.9 from our standard 
run. The dashed line is for binaries with initially circular orbits. 
With the observed stellar radii a collision would occur at about 
r min /a ;£ 0.07. 



Figure 10. Branching ratios as a function of e max . The solid line 
gives the results for encounters in which one binary was ionized. 
The dashed line gives results for a subset of these, namely the 
ones which lead to systems like i Ori (see § 4.1 for the specifics 
of this encounter). 
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Figure 9. Fraction of collisions as a function of e m ax- 
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Figure 11. The fraction of collisions (dotted line) and the 
branching ratio for ionizations (solid curve) as a function of (i. 
Our system of interest is represented with the dashed line. 
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and systems like i Ori which resulted from the standard run. 
There seems to be a slight preference for initial eccentricities 
in the range 0.4 £ e ^ 0.6. 

To further study the dependence of the collision prob- 
ability on the stellar radii we define /3 = d co \\/r, the di- 
mensionless effective radius for collisions. In Figure 11 we 
show the branching ratios for collisions, ionizations and the 
sub-type of t Ori as a function of /3. Changing the stellar 
radii (or 0) by a factor of a few does not effect the collision 
rate much, and therefore the choice of the stellar radii is not 
critical to the results. 



5 THE EVOLUTION OF i ORIONIS AFTER 
THE ENCOUNTER 

Until now we have considered the implications of a binary- 
binary encounter on the dynamical properties of AE Aur, 
/i Col and l Ori. The current observations of the binary 
allow us to constrain its evolution in the ~ 2.5 Myr elapsed 
from the moment of the encounter till now. There are at least 
three further aspects that have to be integrated into a wholly 
plausible model of the interaction under consideration. They 
are: 
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(i) i Ori might be a triple system, and possibly even 
quadruple or quintuple (see § 5.1); 

(ii) stellar evolution models make it difficult to confirm 
that i Ori B is of the same age as any of the other three 
stars (see § 5.2); 

(iii) tidal friction must operate within the present i Ori 
system, and we have to check whether the observed eccen- 
tricity can have been maintained for the 2.5 Myr since the 
collision (see § 5.3). 

We deal with these in turn. 



5.1 Is i Orionis a multiple system? 

The spectroscopic binary i Ori is the brightest (component 
A) of three members of the multiple star ADS 4193. The 
other two members (B and C) are 11" and 50" distant. At 
the distance of t Ori ( ~ 400 pc, sec Table 1) these are quite 
widely separate, and for the time being we discount them as 
gravitationally bound companions. But ADS 4193A is also 
a speckle binary, CHARA 250Aa/Ab, with a separation of 
0.10" (Mason et al. 1998). This suggests an orbital period 
> 40 yr. The brighter component (Aa) of this speckle pair is 
the spectroscopic binary. Strictly speaking we should refer 
to the 09III and B0.8III-IV components of this 29 d binary 
as i Ori Aal and t Ori Aa2, and we do so from now on. 

The close speckle companion raises at least two ques- 
tions. Is it reasonable that the companion was carried along 
during the four-body (or rather 5-body) encounter that we 
have supposed? And could it affect the eccentricity of the 
spectroscopic orbit, by the mechanism of Kozai cycles (Kozai 
1962)? 

We have not attempted to model 5-body encounters 
because of the enormous parameter space that they can be 
drawn from. We can imagine it would be possible for the 
companion to remain bound after an encounter but we doubt 
the energy released would be enough to eject AE Aur and 
H Col with high speed. It might be worth investigating this 
possibility in the future and testing these predictions. 

In a triple system in which the outer orbit is fairly highly 
inclined to the inner orbit, the third body can have a major 
long-term effect on the inner pair, causing its eccentricity to 
cycle between a low and a high value. For instance, if the 
mutual inclination is 60°, e can cycle between and 0.76, or 
between 0.5 and 0.86 (Kiseleva, Eggleton & Mikkola 1998). 
The period of the cycle is ~ Pout /Pin ~ 20000 yr in this case. 
If the triple were formed early in the life of the cluster by 
some random encounter, an inclination of 60°would be quite 
likely; alternatively, the kind of encounter that we think is 
needed to create the binary i Ori Aa would probably ran- 
domize its inclination relative to the (t Ori Aa, i Ori Ab) 
speckle pair. Thus it is possible that the observed high ec- 
centricity of the spectroscopic binary i Ori Aa is not a result 
of the encounter, but of the third body i Ori Ab. 



5.2 Stellar evolution models for AE Aurigae, ji 
Columbae and i Orionis 

We now use stellar evolution model to constrain the ages 
of the four stars. Figure 12 is a theoretical HR diagram 
for masses in the range required. The positions of all four 



Figure 12. Evolutionary tracks in the HR diagram for log 
M=l.l, 1.125,...,1.5. Also shown (thick lines) are isochrones for 
ages 5 Myr and 10 Myr. Sloping nearly-straight dotted lines are 
for log <?=2.8, 3.0,...,3.8, starting from the upper right corner and 
ending near the ZAMS. Approximate positions for l Ori Aal, 
l Ori Aa2, AE Aur and /i Col, based on their temperatures and 
gravities, arc indicated by Aal, Aa2, AE and fj, Col respectively. 
In the scheme the epoch of formation of the eldest binary is la- 
belled as T = 0, the moment of the encounter as T = 4.5 Myr 
and the present epoch as T = 7 Myr. 




stars, according to their effective temperatures and gravi- 
ties as given by Bagnuolo et al. (2001), are indicated. Two 
isochrones are indicated, for 5 and 10 Myr. It is clear that the 
apparent ages of i Ori Aal and t Ori Aa2 are substantially 
different: we estimate 5 and 10 Myr respectively. The uncer- 
tainties appear to be of the order of ±lMyr, but could be 
substantially greater on account of unquantifiable system- 
atic errors. The possibility that discrepant ages can be ac- 
counted for by Roche-lobe overflow can almost certainly be 
discounted by the high eccentricity of the system. The ages 
of AE Aur and n Col are comparable to that of i Ori Aal, 
though AE Aur might be a little younger and [i Col a little 
older. 

The Trapezium Cluster, from which the four stars ap- 
pear to have been ejected, is a well studied young cluster, 
the core of the Orion Nebula Cluster, which contains about 
a dozen OB stars and many more fainter stars. It is esti- 
mated that these stars are g 1 Myr old (Herbig & Terndrup 
1986; Brown, de Geus & de Zeeuw 1994; Prosser et al. 1994). 
Stars later than B2 are often though not always above the 
main sequence. We may therefore wonder if the two giants 
of i Ori are still evolving towards the main sequence rather 
than away from it. This seems improbable, because their 
ejection from the cluster 2.5 Myr ago is amply long enough 
for contraction to the main sequence at their high luminosi- 
ties. This suggests, incidentally, that not all the Trapezium 
stars are as young as 1 Myr, and we appear to require that 
some massive stars have been forming over the last 10 Myr. 
An age spread from 1 to 10 Myr has already been observed 
for stars < 6Mq (Palla & Stahler 1999), and we suggest that 
there may be a similar spread for OB stars. 

An alternative possibility is that the parent cluster is 
one of the older subgroups in Orion such as subgroup la or lc 
(see Brown et al. 1998). These subgroups have a kinematics 
very similar to the one of subgroup Id (which contains the 
Trapezium cluster) and overlap in position on the sky with 
i Ori. 

The uncertainties in the theoretical HR diagram are of 
course substantial, and hard to quantify. The models were 
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computed using the code of Pols et al. (1997), which was 
found to give good agreement with the observational data 
of Andersen (1991). The theoretical uncertainties have to 
be convolved with observational uncertainties but very sub- 
stantial error must be involved if i Ori Aa2 is to be of the 
same age as any of the other three stars. 

A rather long shot is the possibility that the speckle 
companion is in fact the missing lOMyr-old companion of 
i Ori B. This would require a binary-triple interaction of the 
form: 

[i Ori Aa2, i Ori Ab) + (t Ori Aal, (AE Aur, fi Col)) 

((t Ori Aal, u Ori Aa2), i Ori Ab) + AE Aur + p, Col. 

It is hard to assess the likelihood of this but in any event we 
should not ignore the possible existence of a fifth body in a 
complete analysis. 

5.3 The effects of tidal friction on the binary orbit 

An important aspect to consider in the post-encounter evo- 
lution is the effect of tidal friction on the spectroscopic orbit. 
The sum of the radii that we infer is 60% of the periastron 
separation, which would seem sufficient for tidal effects to 
be important. We have used the equilibrium-tide model of 
Hut (1981) as further developed by Eggleton, Kiseleva & 
Hut (1998). This model was successfully tested on the SMC 
radio-pulsar binary 0045-7319 (Eggleton & Kiseleva (2001)). 

The orbital evolution expected since the encounter 
2.5 Myr ago is shown in Figure 13. We find that we could 
reach the present eccentricity and orbital period if we started 
with P=110 days and e=0.91. A high eccentricity is com- 
mon for binaries which undergo exchange-ionization encoun- 
ters as the final eccentricity distribution is expected to be 
thermal (see § 4.2). On the other hand, a long-period binary 
as a final outcome implies an even wider initial target bi- 
nary, which means a lower binding energy available in the 
encounter. It is unlikely to obtain recoil velocities as high 
as 100 km s -1 with a pre-encounter total energy which is 
smaller than the current one. A tight binary (with an or- 
bital period similar to the observed one) is needed in order 
to reproduce the velocities of the single runways. 

Figure 13 shows that parallelization (the alignment of 
the spin axis of the two binary components) and pseudo- 
synchronization took ~ 2 Myr to achieve. Given the uncer- 
tainty in the viscous time-scale - at least a factor of 10 - it is 
possible that they have actually not yet been achieved. The 
pseudo-synchronous rotation period (Hut 1981) is a fraction 
of the orbital period dependent only on eccentricity: the fac- 
tor is 9.51 for e=0.76, i.e. P Iot ~3d currently. With the radii 
computed from the evolutionary tracks of Figure 12, the ro- 
tation speeds we expect are V sini ~ 150 and 115 kms -1 for 
the primary and the secondary respectively. These values 
are somewhat greater than the values (110 and 70 kms -1 ) 
quoted by Marchenko et al. (2000), and may indicate that 
pseudo-synchronism has not yet been reached. 

However we have already noted that the high eccentric- 
ity of the spectroscopic binary may not in fact have been 
generated by the collision but may be an artifact of Kozai 
cycles driven by the third (speckle) body (see § 5.1). In Kozai 
cycles, while the eccentricity fluctuates the semimajor axis is 



unperturbed (to lowest order), and as a result the two com- 
ponents will in the long term seldom be as close together 
as they are at the periastra of the current orbit. Thus tidal 
friction might be negligible in this case. 

The possibility that the eccentricity is due to the third 
body rather than to the interaction could mean that we 
have to revise our earlier conclusion that the apparently 
non-coeval nature of the spectroscopic binary is not due to 
Roche lobe overflow, since the main argument against there 
having been any Roche lobe overflow is that the orbit is 
eccentric. However Kozai cycles are normally suppressed if 
the quadrupolar distortion of the stars in the close binary 
is large; and it usually is large enough as a star approaches 
Roche lobe overflow. Once a star fills its Roche lobe it nor- 
mally grows larger still, until a very late stage of evolution 
when it has lost ~ 70% of its mass. Thus once Kozai cycles 
are suppressed they are unlikely to re-establish themselves 
till this late stage, when the donor would be a small helium 
star rather than a B giant. 

We hope that in the near future more information will 
be available on the speckle orbit, which should move appre- 
ciably in ~10yr. The spectroscopic orbit may also be on 
the margin of detectability by speckle or adaptive optics; a 
direct measure of its inclination could limit the possibilities 
substantially. It will be desirable to search for the speckle 
companion in the spectrum of i Ori. If it can be detected it 
may modify significantly the spectroscopic elements of the 
29 d orbit. 



6 SUMMARY 

We tried to reconstruct the event that ejected the runaways 
AE Aur and /x Col and the binary i Ori from their par- 
ent cluster about 2.5 Myr ago. We here summarize different 
models that could explain the high velocities of the single 
stars and the discrepant ages of the binary components. 

1. The single stars and the binary are all unrelated. As 
shown by Gies and Bolton (1986) and Hoogerwerf et al. 
(2000, 2001), this can be excluded on the basis of their kine- 
matical and evolutionary properties. The same dynamical 
encounter must have ejected the stars from the Trapezium 
about 2.5 Myr ago. In order to explain the high velocities of 
the single stars and the age difference in the binary, the en- 
counter must have involved an ionization and an exchange. 

2. The stars were ejected as a result of a binary-binary 
interaction that occurred in the Trapezium cluster about 
2.5 Myr ago. This model is based on the observation that 
the velocity of the 4-body system is comparable to the ve- 
locity of the Trapezium cluster. The velocities of the stars 
and the period of the binary after the encounter are consis- 
tent with the observed values, so that we can conclude that 
this model can well reproduce the kinematics of the 4-body 
system. If we consider stars with physical radii, a fraction 
of encounters result in a collision between two stars. We 
thus expect to find collision products in dynamically active 
young clusters like the Trapezium. The fraction of collisions 
strongly depends on the initial eccentricity of the two bina- 
ries and therefore we consider it likely that i Ori, AE Aur 
and jx Col were ejected in an encounter between two low- 
eccentricity binaries. In this case, the final properties of the 
system like the recoil velocities of the stars and the orbital 
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Figure 13. The evolution with age of orbital and spin parameters subject to tidal friction. The present epoch is at 2.5 Myr. (a) Left 
panel: eccentricity (decreasing), and the cosine of the inclination £ of the stellar spin axis to the orbital axis (increasing; t Ori Aal only). 
The components were started with rotational axes oblique to the orbital axis; nutation on a short time-scale causes the apparent breadth 
of the inclination curve before parallelization is complete. Interference between the nutation frequency and the rather low data-sampling 
frequency causes some artefacts in the curve, (b) Right panel: orbital period (upper curve), and spin period (lower curve; t Ori Aal only). 
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period and eccentricity of i Ori do not differ from the ones 
obtained for highly eccentric initial binaries. The occurrence 
rate of this type of encounter, once corrected for the collision 
probability, is about one in 10 Myr. 

A possible difficulty in this model is the short circulariza- 
tion time-scale for the orbit of the binary after the encounter. 
Our calculations show that only extremely eccentric post- 
encounter orbits (e ~ 0.9) evolve in an eccentricity like that 
of i Ori after 2.5 Myr. 

In addition, the model relies on the assumption that 
i Ori B and /it Col are in the same binary before the en- 
counter, implying that the two stars must be coeval. Our 
evolutionary calculations show that /i Col might be too hot, 
or equivalently too early in spectral type, to be on the same 
isochrone as t Ori B. However, our knowledge about the tem- 
peratures and gravities of the stars is sufficiently imprecise 
that we can assume coevality for all four stars. 

3. The stars were ejected as a result of a triple-binary in- 
teraction. This model is based on the possibility that t Ori 
is the brightest component of a speckle pair, making it a 
stable hierarchical triple. This would require a 5-body en- 
counter with the same total energy as our standard 4-body 
encounter resulting in the ejection of a triple and two sin- 
gle stars. According to this model, the high eccentricity of 
i Ori may be a result of Kozai cycles driven by the third 
body. We haven't tested this possibility but we can argue 
that this type of encounter would be rather rare and not 
energetic enough to explain the velocities of the runaways. 
It may certainly be worth investigating this scenario to test 
these predictions as well as the ones about corotation and 
pseudo-synchronisation . 



The last word has evidently not been said on the in- 
teresting interaction between i Ori, AE Aur and fi Col, but 
we think it remains a fascinating combination of orbital dy- 
namics, cluster dynamics, stellar evolution and tidal effects. 
In addition to more elaborate numerical simulations, high 
accuracy parameters for the stars are needed in order to 
settle the co-evality issue. New interferometric speckle ob- 
servations, possibly resolving the spectroscopic binary too, 
will be crucial for a better understanding of the dynamical 
encounter that ejected the runaway stars. 
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APPENDIX A: COMPUTATION OF 
CROSS-SECTIONS 

The procedure described in section §3.1 for the computation 
of the maximum impact parameter in each scattering exper- 
iment is incorporated in the actual computation of cross- 
sections. The cross-section for an event X is derived from 
(McMillan & Hut 1996) 



ox 



(Al) 



where Nx,i is the total number of scatterings resulting in 
an outcome X in the shell i and Ni is the total number of 
scatterings performed in the shell i. The uncertainty in the 
cross-section is given by 



2 

i+1 



Ni 



(A2) 



In analogy to the case of binary-single star scattering 
(Hut & Bahcall 1983; Sigurdsson & Phinney 1993) and in 
order to compute physical time-scales for different type of 
interactions, we define a cross-section 



ox = 



ox 



f v oo\' 



(A3) 



7r (a? t + al) 

normalized to the geometric cross-section and corrected for 
gravitational focusing. 

A characteristic time-scale for a process X, defined as 
the mean time between subsequent occurrences of X, can be 
evaluated from (Bacon, Sigurdsson & Davies 1996) 

Tx~ . 2 \ (A4) 
ir (af + a£) V£ n f b o x 

where n represents the mean stellar density in the cluster 
and f b is the fraction of binaries in the core. 



APPENDIX B: CRITERIA FOR STOPPING 
THE INTEGRATION 

During the numerical integration of a binary-binary en- 
counter the status of the system is analysed after each 20 
orbital periods of the initial target binary and the simulation 
is stopped if one of the following criteria is met: 

1. only two stars remain, in which case their orbits are 
determined and the system is classified as a binary if it is 
bound or as two single stars if it is unbound; 
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2. three stars remain, in which case the system is con- 
sidered ionized if the total energy is positive, the stars are 
widely separated and receding from each other while it is 
considered as a binary and a single star if the total energy 
is negative but the single star is escaping from the binary 
center of mass; 

3. the 4-body system is split in binaries whose centers of 
mass are unbound, widely separated and receding. In the 
other cases the system is considered bound and the inte- 
gration continues until a stopping criterion is satisfied or a 
maximum integration time is reached. 



